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ABSTRACT 

We have carried out and analyzed a set of axisymmetric MHD simulations of 
the evolution of a turbulent/diffusive accretion disc around an initially unmagne- 
tized star. The disc is initially threaded by a weak magnetic field where the magnetic 
pressure is significantly less than the kinetic pressure in the disc. The viscosity and 
magnetic diffusivity are modeled by two "alpha" parameters, while the coronal re- 
gion above the disc is treated using ideal MHD. The initial magnetic field is taken 
to consist of three poloidal field loops threading the disc. The motivation for this 
study is to understand the advection of disc matter and magnetic field by the tur- 
bulent/diffusive disc. At early times (< 400 orbits of the inner disc), the innermost 
field loop twists and its field lines become open. The twisting of the opened field 
lines leads to the formation of both an inner collimated, magnetically-dominated jet, 
and at larger distances from the axis a matter dominated uncoUimated wind. For 
later times (> 1000), the strength of the magnetic field decreases owing to field re- 
connection and annihilation in the disc. For the early times, we have derived from 
the simulations both the matter accretion speed in the disc Mm and the accretion 
speed of the magnetic field ub which is determined by measuring the speed of the 
inward motion of the inner 0-point of the magnetic field in the equatorial plane. 
We show that the derived Um agrees approximately with the predictions of a model 
where the accretion speed is the sum of two terms, one due to the disc's viscosity 
(which gives a radial outflow of angular momentum in the disc), and a second due 
to the twisted magnetic field at the disc's surface (which gives a vertical outflow of 
angular momentum). At later times the magnetic contribution to Wm becomes small 
compared to the viscous contribution. For early times we flnd that is larger than 
the magnetic field accretion speed by a factor of ^ 2 for the case where the alpha 
parameters are both equal to 0.1. 

Key words: accretion, accretion discs ~ MHD - black hole physics, magnetic fields, 
jets, stars: winds, outflows 



Early studies of the advection and diffusion of a large- 
scale magnetic field threading a turbulent disc indicated 
that a weak large-scale field would diffuse outward rapidly 
(van Ballegooijen 1989; Lubow, Papaloizou, & Pringie 1994; 
Lovelace, Romanova, & Newman 1994; Lovelace, Newman, 
& Romanova 1997). This rapid outward diffusion may how- 
ever be offset by the highly conducting surface layers of 
the disc where the magnetorotational instability (MRI) and 
associated turbulence is suppressed (Bisnovatyi-Kogan & 
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Lovelace 2007; Rothstein & Lovelace 2008). The magnetic 
field is "frozen-in" in the conducting surface layers which 
tend to fiow inward at approximately the disc accretion 
speed. This conclusion is supported by an analytic model 
for the vertical profiles of the velocity and field components 
of a stationary accretion disc developed by Lovelace, Roth- 
stein, and Bisnovatyi-Kogan (2009). This model predicts 
that the inward or outward transport of the poloidal mag- 
netic flux is determined by both the plasma /3o (the ra- 
tio of the midplane plasma pressure to the midplane mag- 
netic pressure) and the efficiency of the magnetic disc wind 
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in removing angular momentum from the disc (Bisnovatyi- 
Kogan & Lovelace 2012). Guilet and Ogilvie (2012, 2013) 
independently developed an analytic model for the vertical 
structure of a turbulent/diffusive disc threaded by a large 
scale magnetic field, and they find a reduction in the rapid 
outward field diffusion. 

Accretion discs around black holes are considered in 
many cases to be threaded by a large-scale magnetic field 
(Lovelace 1976). This field may be transported inward from 
the interstellar medium by the accreting disc plasma (as in- 
vestigated here), or it may arise from dynamo activity in 
the disc (e.g., Pariev, Colgate, & Finn 2007). Of course the 
discs around magnetized stars may be threaded at large dis- 
tances by the disconnected stellar magnetic field (Lovelace, 
Romanova, & Bisnovatyi-Kogan 1995) . The large-scale field 
may be in the form of magnetic loops threading the disc and 
extending into a low density plasma corona as sketched 
in Figure 1. Differential rotation of the disc acts to open 
the magnetic loops which have footpoints at different radii 
(Newman, Newman, & Lovelace 1992). Also, the differen- 
tial rotation acts to give an axisymmetric field. Such large 
scale magnetic fields can have an essential role in forming 
jets and winds. 

In previous axisymmetric magnetohydrodynamic 
(MHD) simulations of a disc threaded by magnetic loops, 
the disc was treated as a conducting boundary condition 
with plasma outflow and Keplerian azimuthal velocity 
(Romanova et al. 1998). These simulations showed that 
the innermost loop inflates and opens significantly faster 
than the outer loops due to the larger differential rotation 
of the disc close to the star. The opened magnetic fields 
carry away energy, angular momentum, and mass from the 
disc. One or more neutral layers form between the regions 
of oppositely directed magnetic field lines leading to field 
reconnection and annihilation. 

The aim of the present work is to understand the dy- 
namics of the magnetic field loops and the dynamics of the 
disc in response to the twisting of the loops. The magnetic 
field mediates an outflow of energy, angular momentum and 
matter from the disc to a jet and wind. At the same time the 
disc accretion rate can be strongly enhanced by the angu- 
lar momentum outflow to the jet or wind. We treat the disc 
as a viscous/diffusive plasma taking into account fully the 
back reactions of the coronal field on the disc. The turbulent 
viscosity Vt of the disc is modelled with an coefficient 
using the Shakura and Sunyaev (1973) prescription. The 
turbulent magnetic diffusivity r/t is modeled with a second 
a,j coefficient as proposed by Bisnovatyi-Kogan and Ruz- 
maikin (1976). The viscosity and diffusivity are assumed to 
arise from turbulence triggered by the magneto-rotational 
instability inside the disc (Balbus & Hawley 1998), but this 
turbulence is not modeled in the present simulations. The 
low density coronal plasma outside the disc is treated using 
ideal MHD. 

We carry out axisymmetric simulations using a 
Godunov-type scheme to solve the MHD equations, includ- 
ing viscosity and magnetic diffusivity inside the disc as de- 
scribed by Ustyugova et al. (2006). Our initial magnetic 
field configuration consists of three loops in the simulation 
region. New unmagnetized matter is supplied to the disc at 
the outer boundary. 




Figure 1. Sketch of an accretion disc threaded by open and 
closed magnetic field lines from Romanova et al. (1998). Differ- 
ential rotation of the disc will act to rapidly give an axisymmetric 
field configuration considered in this work. 



For this configuration the innermost loop opens up 
rapidly and forms a collimated magnetically dominated jet 
near the z— axis and an uncoUimated matter dominated 
wind at larger distances from the axis. The second loop 
opens at a later time due to the smaller shear in the disc 
and it produces a matter dominated wind. The outermost 
loop reconnects before there is time for it to open. 

The loop configuration allows us to evaluate both the 
accretion speed of the magnetic field and accretion speed of 
the disc matter. This allows a comparison with the analytic 
model of field accretion of Lovelace et al. (2009). 

The paper is organized as follows: We discuss the setup 
for our simulations including the initial and boundary con- 
ditions in Sec. 2. Section 3 describes our results on the dy- 
namics of the disc, the generation of jets and disc winds, 
and the accretion speeds of the matter and magnetic field. 
Section 4 gives the conclusions of this work. 



2 THEORY 

2.1 Basic Equations 

The plasma fiows are assumed to be described by the equa- 
tions of non-relativistic magnetohydrodynamics (MHD). In 
a non-rotating reference frame the equations are 

dp 



dt 
dpw 

dB 

dt 



djpS) 
dt 



+ V • (pv) = , 
+ \/-T = pg , 
c V X E = , 
+ V ■ (pv5) = Q 



(la) 
(lb) 
(Ic) 
(Id) 



Here, p is the mass density, S is the specific entropy, v is the 
flow velocity, B is the magnetic held, T is the momentum 
flux density tensor, E is the electric field, Q is the rate of 
change of entropy per unit volume due to viscous and Ohmic 
heating in the disc, and c is the speed of light. We assume 
that the heating is offset by radiative cooling so that Q = 0. 
Also, g = - [GM/{r - rsf] r (with rg = 2GM/c^) is the 
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gravitational acceleration due to the central mass M with 
the Paczyiiski-Wiita (1980) correction relevant to neutron 
stars. We model the plasma as a non-relativistic ideal gas 
with equation of state 



S = In 



(2) 



where p is the pressure and 7 = 5/3. 

In most of this paper we use spherical (r, 9, 0) coordi- 
nates. However, for some purposes cylindrical coordinates 
are advantageous, and they are denoted [R, cf>, Z). 

Both the viscosity and the magnetic diffusivity of the 
disc plasma are thought to be due to turbulent fluctuations 
of the velocity and magnetic field. Outside of the disc, the 
plasma is considered ideal with negligible viscosity and dif- 
fusivity. The turbulent coefficients are parameterized using 
the a-model of Shakura and Sunyaev (1973). The turbulent 
kinematic viscosity is 



i^t = a. 



(3) 



where Cs is the midplane sound speed, Q,k is the Keple- 
rian angular velocity at the given radii and 5C 1 is a 
dimensionless constant. Similarly, the turbulent magnetic 
diffusivity is 

nt = a^:^ , (4) 

where a,, is another dimensionless constant. The ratio, 

P = - , (5) 

is the magnetic Prandtl number of the turbulence in the disc 
which is expected to be of order unity (Bisnovatyi-Kogan & 
Ruzmaikin 1976). Shearing box simulations of MRI driven 
MHD turbulence in discs indicate that V ~ 1 (Guan & 
Gammie 2009). 

The momentum flux density tensor is given by 



Tik = pSik + pviVk + 



8^' 



Bi Bk 



(6) 



where Tik is the viscous stress contribution from the tur- 
bulent fluctuations of the velocity and magnetic field. As 
mentioned we assume that these can be represented in the 
same way as the collisional viscosity by substitution of the 
turbulent viscosity. Moreover, we assume that the viscous 
stress is determined mainly by the gradient of the angu- 
lar velocity because the azimuthal velocity is the dominant 
velocity of the disc. The leading order contribution to the 
momentum fiux density from turbulence is therefore 



Tr4, — —Vtprsvat 



T04, 



dr 



duj 



(7a) 



(7b) 



where lj — v^i/r sin 9 is the plasma angular velocity. 

The transition from the viscous-diffusive disc to the 
ideal plasma corona is handled by multiplying the viscosity 
and diffusivity by a dimensionless factor ^(p) which varies 
fi-om ^ = for p s; 0.25pd to ^ = (4/3) (p/pd - 0.25) for 
0.25pti < p < pd to ^ = 1 for p > pd (see Appendix B of Lii, 
Romanova, & Lovelace 2012). The disc half-thickness h is 



taken to be the vertical distance from Z = to the 0.5pd 
surface. 



2.2 Initial Conditions 

2.2.1 Initial Magnetic Field 

The initial magnetic field is described in cylindrical coordi- 
nates. This field is taken to be force-free in the sense that 
J X B = in the region \Z\ > with B = {Br,B^,Bz), 
where Br = -R-^d^/dZ and Bz = R'^d^/dR. Here, 
'i!{R, Z) is the fiux function which labels the field lines, 
B ■ V'l' = 0. This function satisfies the Grad-Shafranov 
equation. 



A**(i?,Z) = -H(*) 



(8) 



where 



A* 



dR^ 



RdR'^ dZ^ 



and H = H{^) = RB^(R,Z) is the poloidal current 
function (Lovelace et al. 1986). For simplicity we take the 
poloidal current to be proportional to the fiux, H{^) = , 
where fc is a constant (Newman et al. 1992). The relevant 
solution to equation (8) is 



■^{R,Z) = AiRJi{aiR)e-^'^^^ +A2RJi{a2R)e' 



, (9) 



where Ji is a Bessel function of the first kind, Ai, A2 are in- 
tegration constants and ai = (k^ +6?)^^'^ for i = 1, 2. Qual- 
itatively, the field appears as a number of loops threading 
the accretion disc. We have chosen the parameters ai such 
that three loops fit in our simulation region. 

The solution (9) is valid in the region IZj > and 
assumes that initially there is a thin current carrying disc 
in the Z — plane. The cusp in the initial field at Z = 
disappears rapidly in a time tinit ~ h^/Vt = {an^K)~^ = 
Po(27rQ:^)~^(r/ro)^''^ due to the diffusivity of the disc. Here, 
ro is the reference radius and Pq period of the Keplerian 
orbit at this radius as discussed in §2.4. We have assumed 
h/R ^ Cs/vk which neglects the magnetic compression of 
the disc discussed by Wang, Sulkanen, and Lovelace (1999). 
For most of the range of r/ro this time is much smaller than 
the field evolution time scale. 



2.2.2 Matter Distribution 

Initially the matter of the disc and corona are assumed to 
be in mechanical equilibrium (Romanova et al. 2002). The 
initial density distribution is taken to be barotropic with 

^^^^ _ fp/Tdisc p>Pb and rshi9^n, ^^^^ 
\p/Tcoi P < Pb or rsin6'^r-6, 

where pb is the level surface of pressure that separates the 
cold matter of the disc from the hot matter of the corona 
and rt is the initial value of the inner radius of the disc. 
At this surface the density has an initial step discontinuity 
from value p/Tdisc to p/T^ov. 

Because the density distribution is barotropic, the ini- 
tial angular velocity is a constant on coaxial cylindrical sur- 
faces about the axis. Consequently, the pressure can be 
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determined from Bernoulli's equation, 

f (p) + "I> + $c = const , (11) 

where <& = —GM/\r — rc\ is the gravitational potential with 
the Paczyhski-Wiita correction, $c = X-^ingC'^C ^^{0 is 
the centrifugal potential, which depends only on cylindrical 
radius R = rsin^, and 

p(^p^ _ Iodise ln(p/Pb) P> Pb and rsinO^n , ^^^-j 
[Tcor ln(p/pt,) p < Pb or rsin6'^r-6. 

^.2.5 Angular Velocity 

Initially the inner edge of the disc is located at ri, = 5ro in 
the equatorial plane, where ro is a reference length discussed 
below. The initial angular velocity of the disc is slightly sub- 
Keplerian, 

n|e=v2 = (1 - 0.003)0^(7-) r>rb, (13) 

Inside of ri,, the matter rotates rigidly with angular velocity 

= (l- 0.003)t7A'(r-i,) r sC r^. (14) 

The corotation radius r^r is the radius where the angular 
velocity of the disc equals that of the star; that is, Vcr = 
(GMt . In this study we have chosen this radius to 
be the initial inner radius of the disc with Va = 5ro. 

2.3 Boundary Conditions 

Our simulation region has four boundaries: The surface of 
the star, the midplane of the disc, the rotation axis, and the 
external boundary. For each dynamical variable we impose 
boundary conditions consistent with our physical assump- 
tions. 

We assume axisymmetry as well as symmetry about 
the equatorial plane. On the star and the external bound- 
ary we want to allow fluxes and so impose free boundary 
conditions dT jdr = 0, where F are the dynamical vari- 
ables. In addition, along the external boundary in the disc 
region 9 = 72° — 90°, we allow matter to inflow but the in- 
flowing matter has zero magnetic flux. In the coronal region 
= 0° — 72° we allow matter, entropy and magnetic flux to 
exit the simulation region. 

In addition to the boundary conditions, which are re- 
quired for the well-posedness of our problem, we impose 
additional conditions on variables to eliminate numerical 
artifacts in the simulations. For instance, we require that 
the radial velocity at the surface of the star be negative. 
Thus there there is no outflow of matter, angular momen- 
tum, or energy from the star. We also require that along 
the external boundary the radial velocity be inwards inside 
the disc and outwards outside the disc so the disc matter 
will tend to accrete and matter in the corona will tend to 
be ejected. 

Figure 2 shows the grid used in the described simu- 
lations. There are Ne = 30 constant width cells in the 
0— direction. In the r— direction, the A''^ = 67 grid cells in- 
crease in width as drj+i = (1-1-0. 0523)drj so as to give curv- 
linear rectangles with approximately equal sides (Ustyugova 
et al. 2006). The dependence of our results on the grid res- 
olution is discussed in Appendix A. 




Figure 2. Grid used in the simulations. 



Parameters 


Symbol 


Value 


mass 


M, 


2.8 X 10^^ g 


length 


ro 


1.0 X lO^cm 


magnetic field 


Bo 


10* G 


time 


Po 


4.6 X 10-*s 


velocity 


Vo 


1.4 X lOiOcm/s 


density 


PO 


1.3 X 10-5 g/cm3 


accretion rate 


M. 


2.8 X lO"'' M0/yr 


disc power 


Eo 


1.66 X 10^'^ erg/s 



Table 1. Mass, length, and magnetic field scales of interest and 
the corresponding scales of other derived quantities 



2.4 Dimensional Variables 

The MHD equations are written in dimensionless form so 
that the simulation results can be applied to different types 
of stars. The mass of the central star is taken as the ref- 
erence unit of mass. Mo = M«. The reference length, ro, 
is taken to be half the radius of the star. The initial in- 
ner radius of the disc is rj, = 5ro. The reference value 
for the velocity is the Keplerian velocity at the radius ro, 
Vo = (GMo/ro)^^'^ ■ The dimensionless temperature is T/vq. 
The reference time-scale is the period of rotation at ro, 
Po ~ 2-Kro/vo. From the MHD equations, we get the re- 
lation poVo = Bo, where Bo is a reference magnetic field 
and Po is a reference density both at ro. We take the refer- 
ence magnetic field Bo to be such that the reference density 
is appropriate for the considered star. The reference mass 
accretion rate is Mo = povoro- The reference disc accretion 
power is Eo — GMoMo/ro- The initial dimensionless tem- 
perature in the disc is Tdisc ~ (p/p)disc = 5 x 10"'', and the 
initial temperature in the corona is Tcor ~ {p/p)coi = 0.5. 

Results obtained in dimensionless form can be applied 
to objects with widely different sizes and masses. However, 
the present work focuses on neutron stars with the typical 
values shown in Table 1. 
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Figure 3. Plot of poloidal magnetic field lines (white) and matter density (colour) for oiu = 0.1 and a,, = 0.1 at t = 0, = 300, .., 1500, 
where t is the time measured in units of Pq which is the period of the Keplerian orbit at the reference radius ro (see §2.4). By t = 300 
sufficient matter has fallen into the star to drag in the magnetic field. The latter forms a well coUimated jet along the axis and a wind 
alone thcLdigp. .Ttvft Pff fi?)')^ fo n jo me time, but eventually decays due to magnetic field reconnection and annihilation. 
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3 RESULTS 

We have carried out a large number of simulation runs for 
different values of the viscosity 0.05 ^ ^ 0.3 and diffu- 
sivity 0.01 ^ an ^ 0.3 parameters and find that the simula- 
tions exhibit similar qualitative behaviour. Figure 3 shows 
the evolution of the poloidal field projections for a repre- 
sentative case where =0.1 and = 0.1. The field lines 
of the innermost loop are pulled in towards the star by the 
accreting disc matter. When this loop reaches the star's 
surface it opens up. The inner half of the loop extends ver- 
tically upwards from the star, supporting a magnetically 
dominated jet along the z— axis. The outer half of the loop 
threading the disc projects outwards from the disc at about 



45° to the disc normal, and it supports a magnetic disc 
wind. The middle and outer magnetic loops move inward 
only gradually, and they decrease in strength due to field 
line annihilation inside the disc. 

Figure 4 shows the radial dependence of the inverse 
"plasma beta" which is the ratio of the magnetic pressure 
to the plasma pressure at the disc midplane, 



8ttp{R, 0) 



(15) 



where Bzo = Bz{R, Z = 0). Note that /Jg"^ = {vAo/csof/2, 
where vao = Bzo/ ^^np is the midplane Alfven speed and 
Cso is the midplane sound speed. For the assumed symme- 
try about the equatorial plane Bz is the only non-vanishing 




Figure 4. The ratio of magnetic pressure to the plasma pressure 
fi^^ in the disc midplane at t = (full), t = 100 (dashed) and 
t = 300 (dot-dashed) for ai, = 0.1 and a,, = 0.1. 




© 2012 RAS, MNRAS 000,[l]{l4] 



Advection of Matter and B-Fields in Alpha-Discs 7 



field component at Z = 0. Initially 13^^ is significantly less 
than unity over most of the disc [R > 5). Consequently, 
the magneto-rotational instability would be expected to oc- 
cur in an actual disc with the same parameters (Balbus & 
Hawley 1998). At later times, [Bz{R,0)]^ increases but the 
plasma pressure also increases so that the change in 
are not very large. 

3.1 Opening of Field Loops and Field 
Annihilation 

The initial three poloidal field loops {t — panel of Fig- 
ure 3), have footpoints at the approximate radii Rk = 
5, 20, 35, & 50 (fc = 1,..,4). The time-dependent open- 



ing of large-scale magnetic field loops with footpoints at 
different radii in a Keplerian disc has been well analyzed 
(Newman, Newman, & Lovelace 1992; Lynden-Bell & Boily 
1994; Romanova et al. 1998; Ustyugova et al. 2000; Lovelace 
et al. 2002). We use the numerical condition of Lynden- 
Bell and Boily (1994) that a field loop opens after there 
is a differential rotation of the footpoints by > 3.63 radi- 
ans. For each loop, the differential rotation in a time is 
A(p = tkl^KiRk) - ^^(-Rfc+i)]. With = 3.63 we obtain 
the opening times for the three loops, t\ = 40.8, ti — 329, 
and = 760 in our dimensionless units. The observed rapid 
opening of the first loop (Figure 3) agrees qualitatively with 
ti. The fact the second loop does not open in a time ~ t2 
may be explained by the overlying magnetic field of the first 




Figure 6. The magnetic field Bz{R, Z = 0) along the disc mid- 
plane at t = (full), t = 200 (dashed) and t = 400 (dot-dashed) 
for Ov = 0.1 and a,, = 0.1. 
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I I I I I I 
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10 20 30 40 50 
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Figure 7. The toroidal magnetic field B^{R, Z = h) at the disc 
surface at t = (full), t = 200 (dashed) and t = 400 (dot-dashed) 
for a,y = 0.1 and = 0.1. 
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loop. The long time ts required for the opening of the third 
loop means the diffusion of the magnetic field has sufficient 
time to cause significant field annihilation. 



The time scale for the field to diffuse over a dis- 
tance ARk/2 — {Rk+i — Rk)/'2 can be estimated as ~ 
(A7ife)^/(4r;t), where jyt is evaluated at Rk — {Rk+Rk+i)/^. 
For ar, = 0.1, we find n « 2500, T2 ^ 1700, and ts ^ 1400. 
Note that rg for the outer loop is less than the duration 
of our runs so that the magnetic field decays significantly 
before the end of the runs. 



3.2 Matter Advection in the Disc 

Figure 5 shows the average accretion speed of the disc 
matter as a function of R at different times, where 

Mm = dZp(R,Z)vR{R,Z) . (16) 

(y(R) J-h 



Here, 



cj{R) = / dZp{R, Z) , 

-h 



(17) 



is the surface mass density of the disc. The conservation of 
mass gives 

9(i?Cr) d{RGUm) _ 1 5M^ 



dt 



dR 



TT SR ' 



(18) 




Figure 8. The radial component of the magnetic field 
Bii{R, Z = h) at the disc surface at t = (full), t = 200 (dashed) 
and t = 400 (dot-dashed) for = 0.1 and a,, = 0.1. This quan- 
tity has an important role in determining the diffusive advection 
of the magnetic field Jig^ (equation 26). 
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Matter advection in the disc is measured by Um{R,t), 
which is determined by the disc's turbulent viscosity (which 
causes the radial outflow of angular momentum) and by 
the torque of the large-scale magnetic field (which causes a 
vertical outflow of angular momentum to magnetic jets or 
winds). This is described by a simple analytic model where 



Mm = Umv + UmB , (19) 




Figure 10. Three dimensional view of a magnetic field line orig- 
inating from the disc at ij = 6 at time t = 300. The twist of the 
field line is such that the field transports angular momentum out 
of the disc. That is, BzB^ <0. 




Figure 11. The solid curve shows the average accretion speed 
of the disc matter at the maximum of Bz{R,Z = 0) from the 
simulations for = 0.1 = a,,. The initial radius of the maxi- 
mum is R = 22. The dashed curve is from the advection model 
(eqns. 19 and 20)) with parameters C\ = 0.05 and C2 = 0.2. 
The model breaks down when the maximum approaches the star 
and interacts strongly with it. This happens at different times 
for different viscosity and diffusivity values. 



where SM^/SR = 2irR{pvz)z^h is the mass outflow rate 
from the top surface of the disc to the wind. The mass 
accretion rate to the star from the upper half-space is M* = 
(7r7icrUm)fl=ro ■ We find that the ratio of the time-averaged 
mass loss rate of the wind is typically small compared to 
M*. 

Figure 6 shows the midplane magnetic field Bz{R,0) 
at different times. Because of the assumed symmetry of the 
magnetic field about the equatorial plane, this is the only 
non-zero field component at Z — 0. 
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400 



Figure 12. Accretion speed of the magnetic field ug in the 
equatorial plane derived from the motion of the innermost zero 
crossing of Bz {R, Z = 0,t) and the matter accretion speed Um 
at the same location for = 0.1 = a,,. The initial radius of 
zero crossing is Ro{t = 0) = 13.7. 



where 



VK, UmB 



-C'zB^hBzo 



(20) 



(Lovelace et al. 2009; 1994), where h is the half-thickness 
of the disc, vk = {GM / R)^^^ is the Keplerian velocity, a 
is the surface mass density of the disc, B^h ~ B^{R,h) 
is the toroidal magnetic field at the disc surface, Bzo = 
Bz{R,Z = 0) is the midplane magnetic field, and Ci and 
C2 are dimensionless constants of the order of unity. The 
first term represents the accretion speed contribution due 
to the turbulent viscosity of the disc, while the second term 
the accretion speed contribution due to the outflow of an- 
gular momentum from the disc surfaces due to the twisted 
magnetic field in the corona. 

Figure 7 shows the profiles of the toroidal magnetic 
field at the disc surface Btf,{R,h) at a sequence of times. 
This quantity is important for the outfiow of angular mo- 
mentum from the disc which in turn determines the mag- 
netic contribution to the accretion speed UmB in equation 
(20). 

Figure 8 shows the radial profiles of Bjih at a sequence 
of times. This quantity is important for the radial diffusion 
of the magnetic field as discussed below in §3.3 Figure 9 
shows the radial variation of h/ R at a sequence of times. 

Figure 10 shows that the twist of a sample magnetic 
field line above the disc is such that B^hBzo < which 
corresponds the to outflow of angular momentum from the 
disc. 

We flnd that during early times (i < 400), the Bz{R, 0) 
field "propagates" inwards towards the star as shown in 
Figure 6. We study this motion by measuring the positions 
of the inner maximum of Bz{R,0) as a function of time. 
The maximum moves inwards until interactions with the 
star cause a more complicated behavior. At the Bz{R,0) 
maximum at a given time t, we calculate the accretion speed 
UiniR,t) in our simulation data using equation (16). We 



also calculate B^h{R,t), Bzo{R,t), h/R, and a{R,t) which 
permits us to compare the observed accretion speed with 
the prediction of the advection model (equations 19 and 
20). We find that the reasonable values of Ci — 0.4 and 
C2 — 0.2 give satisfactory agreement between the model 
and our various simulation runs for t < 400. 

Figure 11 shows the model and the measured simula- 
tion accretion speeds for sample cases. 

3.3 Magnetic Field Advection in the Disc 

The advection of the poloidal magnetic field is described by 
the equation 



d{RBz) d{RBzUBi) 



dt 



dR 



d_ 
dR 



ritRBRh 



+ VtR 



dBz 
dR 



(21) 

(Lovelace et al. 2004), where Bru = Br{R,Z = h). For 
simplicity we have neglected terms of order \dh/dR\ relative 
to unity. Here, 



UBi 



dZBz{R,Z)vR(R,Z) 



dZBz{R, Z), 



(22) 

is the magnetic field advection speed of an ideal, perfectly 
conducting disc (774 = 0). Note that the matter advection 
speed density weighted average over the disc thick- 

ness of —vr whereas UBi is an average of —vr weighted by 
Bz- For smooth profiles of p and Bz the two speeds will be 
be comparable. 

The vertical magnetic flux threading the disc inside the 
first 0-point where _Bz(_Ro,0) — (or between successive 
0-points) decreases in general with time due to the diffu- 
sivity. From equation (21) we have 

fflo 



d 
dt 



RdRBz 



VtRBE 
h 



+ r)tR 



dBz 
dR 



(23) 



R^Ro 

For example, for Bz < inside Ro, both terms on the 
right-hand side of equation (23) are seen to be positive so 
that the magnitude of the fiux decreases. 

For rjt > the field advection speed is the sum of the 
ideal and diffusive contributions, 



with 



Here 



UB = UBi + UB-n 



d{RBz) d{RBzUB) 



dt 



UB-n = 



dR 



= 



ntBRh ^ 77* dBz 



hBi 



Bz dR 



(24) 



(25) 



(26) 



is the diffusive advection speed. 

From our simulation data we can calculate the advec- 
tion speed of the magnetic field UBi by tracking the location 
Ro{t) of the Bz{R,Z — 0) zero crossings which occur at 
"O-points" of the poloidal magnetic field Bp. Close to the 
zero crossing, Bz{R,0) = const(_R — Ro) is an odd func- 
tion of R — Rq. Furthermore, Bz is also an odd function 
of i? — Ro ■ Thus UBrj is an odd function about the 0-point 
proportional to {R — Ro)~^ because of the Bz denomina- 
tors in equation (26). The magnetic field moves symmetri- 
cally inward towards the 0-point where it annihilates. The 



© 2012 RAS, MNRAS 000,[l]{l4] 



Advection of Matter and B-Fields in Alpha-Discs 11 



mathematical singularity is smoothed out by the finite grid. 
The magnetic field inside a current-carrying resistive wire 
behaves in the same way. Consequently, at J? = Ro we have 
wsrj ~ 0. The diffusivity has no influence on the motion of 
the O-point. That is, dRo/dt — —ub = —UBi- 

We can compare field advection speed ub (at an O- 
point) to the matter advection speed Um at the same lo- 
cation. As mentioned these two speeds are expected to be 
comparable. 

Figure 12 shows sample comparisons oi ub and iim. 
For a smaller difi'usivity with fixed viscosity the dilTerence 
between matter and field accretion speeds is smaller. For 
a larger viscosity relative to diffusivity, Um is significanty 
larger than ub- 



3.4 Late times 

At late times {t > 1000), the magnetic field decays appre- 
ciably owing to reconnection and field annihilation. Conse- 
quently the accretion speed is due mainly to the disc vis- 
cosity. 

The mass accretion rate to the star from the top half space 
is M* = 7r(JiEuni)*, where E is the surface mass density 
of the disc and the asterisk subscript indicates evaluation 
outside the star. 

Figure 13 shows the time dependence of the accretion 
rate to the star and the mass outflow rate in the wind for 
two viscosity values and = 0.1. For 0.05 ^ ^ 0.3 
we find that Af« at late times (t > 1000) is approximately 
proportional to At late times M» is independent of the 
diffusivity a^. 



3.5 Jet and Wind 

In this work we observe both a collimated jet along the 
Z— axis and an uncollimated disc wind. 



3.5.1 Jet 

The fluxes of angular momentum and energy through the 
spherical surface [r = 44, 0° ^ ^ 21°] are shown in Figure 
14 The angular momentum flux can be separated into a part 
from the matter and a part due to the magnetic field, 



L — Lm + Lf 



dS ■ I pr sin{9)v^'Vp 



rsm{e)B^Bf 



(28) 

Similarly, the energy flux can be separated into a contribu- 
tions carried by the matter and that carried by the Poynting 
flux, 



E — Em + Ef 



dS 



ExB 



(29) 



The jet is strongly dominated by the electromagnetic field: 
The angular momentum flux is carried predominantly by 
the magnetic field and the energy flux is carried predom- 
inantly by the Poynting flux. Such jets were hypothesized 
by Lovelace (1976) and first observed in axisymmetric MHD 
simulations by Ustyugova et al. (2000). 



3.5.2 Disc Wind 

The rates of energy, angular momentum and mass flux 
through the surface [r = 44, 21° < < 72°] is shown 
in Figure 15. We have chosen the upper bound for 6 by 
requiring that the wind stay outside the disc. 

Figure 14 shows the different components of the wind 
angular momentum flux and the components of the energy 
flux. The angular momentum and energy fluxes are dom- 
inated by the matter components. This is the opposite of 
the case for the jet. 

Figure 14 shows the jet and wind total energy fluxes 
normalized to the "accretion power" -Eacc = G'M«M«/(2rt). 
The large initial values of the ratios results from the fact 
that M* is initially zero. Note that the two ratios are com- 
parable. 



4 CONCLUSIONS 

We have analyzed a set of axisymmetric MHD simulations 
of the evolution of a turbulent/diffusive accretion disc ini- 
tially threaded by a weak magnetic field with midplane 
plasma beta /3o is significantly larger than unity. The vis- 
cosity and magnetic diffusivity are modeled by two a pa- 
rameters, one for the viscosity and the other for the 
diffusivity a,,. The coronal region above the disc is treated 
using ideal MHD. The initial magnetic field is taken to con- 
sist of three poloidal field loops threading the disc between 
its initial inner radius and to its ten times larger outer ra- 
dius. This field configuration allows the derivation of the 
advection speed of the magnetic field. 

Recent theoretical studies discussed the importance 
of the magnetic field extending from a turbulent disc 
into a low density non-turbulent /highly conducting corona 
(Bisnovatyi-Kogan & Lovelace 2007; Rothstein & Lovelace 
2008; Lovelace et al. 2009; Bisnovatyi-Kogan & Lovelace 
2012; Guilet & Ogilvie 2012, 2013). These treatments aU 
considered stationary or quasi-stationary conditions and a 
disc threaded by a poloidal magnetic field of a single polar- 
ity. In contrast the simulations discussed here are strongly 
time-dependent and involve multiple poloidal field polari- 
ties in different regions of the disc. Consequently, a direct 
comparison of the theory and simulations is not possible. 
The simulations clearly show the inward advection of the 
magnetic field at about the same speed as the matter ad- 
vection before the field decays by annihilation. 

At early times {t < 400), we find that the innermost 
field loop twists and its field lines become open. For the 
different field loops we estimate two important time scales: 
One is the time scale for each loop to open due to differ- 
ential rotation of its foot points, and the other is the field 
annihilation time scale owing to the disc's magnetic diffu- 
sivity. The innermost field loop opens rapidly before there 
is significant annihilation. On the other hand the outer loop 
decays significantly before there is time for it to open. The 
twisting of the opened field lines of the inner loop leads 
to the formation of both an inner collimated magnetically 
dominated jet and at larger distances from the axis a matter 
dominated uncollimated wind. For later times (> 1000), the 
strength of the magnetic field decreases owing to field re- 
connection and annihilation in the disc. For the early times, 
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Figure 14. - Left panel: The jet angular momentum flux carried by the matter (dashed curve), the magnetic field (dot-dashed curve), 
and the total flux (solid curve) for = 0.1 and = 0.1. Right-panel: The jet energy flux carried by the matter (dashed curve), the 
Poynting flux (dot-dashed curve), and the total flux (solid curve) for Oi/ = 0.1 and q,, = 0.1. The jet is strongly dominated by the 
Poynting flux. 



we have derived from the simulations both the matter ac- 
cretion speed in the disc and the accretion speed of 
the magnetic field ub- We show that the derived u-m agrees 
approximately with the predictions of a model where the 
accretion speed is the sum of a contribution due the disc's 
viscosity (which gives a radial outflow of angular momen- 
tum in the disc) and a term due to the twisted magnetic 
field at the disc's surface (which gives a vertical outflow of 
angular momentum) (Lovelace et al. 2009; 1994). At later 
times the magnetic contribution to Um becomes small com- 
pared with the viscous contribution. Also for early times 
we find that is larger than the magnetic field accretion 
speed it_g by a factor ~ 2 for the case where = 0.1 = a,,. 
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and the total (solid curve) for = 0.1 and On = 0.1. Right panel: The vfind energy flux carried by the matter (dashed curve), the 
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ari = 0.1. Right panel: The total wind power also normalized by the accretion power Eacc for 
values of the ratios is due to the fact that the M« is initially zero, 
partition of the range of 6 into wind and disc regions. 



GM,M,/{2r,) for = 0.1 and 



: 0.1 and 



: 0.1. The large initial 



The fact that -Ewind is temporarily negative is due to the fixed 



Gullet, J. & Ogilvie, 
| |arXiv :1212.0855) 

M.M., 



G.I. 2013, MNRAS, 



press 



Lii, P., Romanova, 

420, 2020 
Lovelace, R.V.E. 1976, Nature 
Lovelace, R. V. E., Mehanian, 

M.E. 1986, ApJS, 62, 1 
Lovelace, R.V.E. , Li, H., Koldoba, A.V. 

manova, M.M. 2002, ApJ, 572, 445 
Lovelace, R.V.E., Romanova, M.M., & 

1995, MNRAS, 275, 244 
Lovelace, R. V. E., Newman, W. I 

ApJ, 484, 628 
Lovelace, R.V.E., Romanova, M.M. 
437, 136 

Lovelace, R.V.E., Rothstein, D.M 
2009, ApJ, 701, 885 



& Lovelace, R.V.E. 2012, MNRAS, 



262, 649 

C., Mobarry, C.M., 



Sulkanen, 



: Ro- 



Ustyugova, G.S., 
Bisnovatyi-Kogan, 
& Romanova, M. M. 
k Newman, W.I. 1994, ApJ 
& Bisnovatyi-Kogan, G.S. 



G.S. 



1997, 



Lubow, S. H., Papaloizou, J. C. B., & Pringle, J. E. 1994, MN- 
RAS, 267, 235 
Lynden-BeU, D., & Boily, C. 1994, MNRAS, 267, 146 
Newman, W.I., Newman, A.L., & Lovelace, R.V.E. 1992, ApJ, 
392, 622 

Paczynski, B., & Wiita, P. 1980, A&A, 88, 23 

Pariev, V.I., Colgate, S.A., Finn J.M. 2007, ApJ 658: 129-160. 

Romanova, M.M, Ustyugova, G.V., Koldoba, A.V., Chechetkin, 

V.M., & Lovelace, R.V.E., 1998, ApJ 500, 703 
Romanova, M.M, Ustyugova, G.V., Koldoba, A.V., Chechetkin, 

V.M., feLovelace, R.V.E. 2002, ApJ, 578, 420 
Rothstein, D.M., & Lovelace, R.V.E. 2008, ApJ, 677, 1221 
Shakura, N.I., & Sunyaev, R.A. 1973, A&A, 24, 337 
Ustyugova, G.V., Lovelace, R.V.E., Romanova, M.M., Li, H., & 

Colgate, S.A. 2000, ApJ, 541, L21 
Ustyugova, G.V., Koldoba, A.V., Romanova, M.M, & Lovelace, 

R.V.E., 2006, ApJ 646:304-318. 



© 2012 RAS, MNRAS 000,[T|fT4l 



14 S. Dyda et al. 




10 20 30 40 50 
R 

Figure Al. Radial dependences of Bz{R, Z = 0) at a sequence 
of times for the higher (soUd curves) and lower resolution runs. 
For both cases a„ = 0.1 and a,, = 0.1. 
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APPENDIX A: DEPENDENCE ON GRID 
RESOLUTION 

We have tested the dependence of our results on the grid 
by running higher resolution cases compared with resolu- 
tion used for this study, (A^e, Ng =) = (31, 67). We 
have run cases with (41, 87) and (50, 100). Figure Al 
shows the radial dependence of the mid plane magnetic field 
Bz{R, Z = 0) at t = 300 for the three grid resolutions. The 
radius of the first zero crossing of Bz{R,0) decreases by 
about 14% going from the low to the intermediate reso- 
lution. It decreases by a further 2% going from the inter- 
mediate to the high resolution grid. Thus the convergence 
is rapid. This indicates the accuracy of the field advection 
speed Ms in Figure 12 and suggests that the actual speed 
is higher by about 16%. 
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